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ABSTRACT 


The problem of orbital parameter estimation using angles only observations is 
examined. Direction cosine measurements, obtained from satellite passage of an 
earth-based stationary planar radar beam, are assimilated by an extended Kalman 
filter to improve estimates of a classical orbital element set. Several progressively 
comprehensive orbital motion models are considered and compared. The relative 
effectiveness of these models is illustrated by applying them to actual satellite 


data. 
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Il. INTRODUCTION 


Earth-orbiting satellites are catalogued in the form of orbital element sets. For 
maximum usefulness of these catalogs to military and civil space systems, element 
sets need to be updated as frequently and as accurately as possible. A nonlinear 
filter was designed to perform these updates in near-real-time using observations 
of satellites piercing a planar radar beam. 

The filter predicts time-varying orbital parameters according to a non-linear 
orbital dynamics model. This model includes, along with two-body Keplerian 
motion, first order oblateness effects, both secular and periodic. simple decreasing 
exponential atmospheric drag effects, and white noise on both the dynamic model 
and the measured observations. This simple model is found sufficient to track 
several selected satellites, with proper tuning of the filter. It should be noted that 
even this simple non-linear model gives rise to comparatively complex filter 
expressions. 

Applying this filter to the orbit improvement problem helps meet the need for 
continual determination of orbit decay, collision avoidance, satellite maneuvers, 


etc. Current systems use a daily batch least squares differential correction method. 


Il. REVIEW OF LITERATURE 


This research involves the integration of a simple, yet adequate, orbital 
dynamics model with a necessarily suboptimal extended Kalman estimator. This 
combination is applied to a geometrically unique detection scheme to provide near 
real-time improvement of orbital elements for selected satellites. The accuracy 
level under consideration falls under the umbrella of the general perturbations 
problem of orbital mechanics, which 1s covered by numerous textbooks, a 
sampling of which are referred to here [1-6]. In general, solution procedures 
involve a gravitational potential function (earth oblateness, third body, etc.) 
describing its contribution to the perturbation of a classical two-body Keplerian 
dynamic model. The resultant expressions take the form of variations in some set 
of orbital parameters. The two broad divisions of orbital parameter representation 
in the literature are: 1) some set of 6 orbital elements which describe the orbital 
path as a conic section and its orientation in inertial space, and 2) Cartesian 
position and velocity of the orbiting body. The 6 classical orbital elements were 
chosen for their physical significance in terms of how the satellite's orbit is being 
affected by various perturbing forces. The particular satellite propagation theory 
being used here is loosely based on Brouwer's theory [7-11] for a future 
comparison basis with existing orbit improvement procedures. The detection 
scheme 1s an earth-fixed planar radar, the crossing of which causes power to be 


reflected back to one or more of multiple earth-based receivers [12]. 


The most poorly modeled phenomenon affecting low earth orbits is the 
problem of atmospheric drag. The overwhelming consensus [13-15] is that even 
the most complex atmospheric density models can only be consistently accurate to 
within about 20%. The Kalman filter incorporates a very simple decreasing 
exponential density model [16] corrupted by noise to track drag perturbations. 

The extended Kalman filter is used in many situations for nonlinear parameter 
estimation [17-19] and has been applied to the field of orbit estimation in the form 
of orbital parameter improvement by differential correction techniques [20,21]. It 
is known that such a filter works well when observational data are plentiful 
[22,23]. However, proper application of such a filter to a highly nonlinear, 
minimal observation scenano, such as that posed by the fixed radar plane detection 


scheme, can also successfully "track" satellites. 


III. TWO-BODY ORBIT FORMULATION 


Two-body orbital motion for a spherical earth can be described by a set of six 
parameters (three-dimensional, second order problem). Traditionally, these 
parameters take one of two forms. One method describes satellite motion in terms 
of Cartesian position and velocitv. These vectors are propagated according to the 
two-body non-linear equations of motion. The other takes advantage of the fact 
that the equations of motion describe an orbital path which is a conic section, 
which, for earth-orbiting satellites, is always an ellipse. Satellite motion can then 
be analyzed as a mean angular displacement along that elliptical path. The second 
method is used here because of the linear nature of the mean angular motion and 


its claritv of phvsical significance. 


A. KEPLER TWO-BODY DYNAMICS 

The orbital path is described as an ellipse of fixed shape and size located in a 
plane of fixed orientation with respect to the earth's inertial reference frame. The 
quantities describing this path, the classical orbital elements [1], consist of five 
constants and one time-dependent quantity. Two of the constants (a, e) describe 
size and shape of the orbit, while the other three (1, Q, œ) orient the orbital plane in 
space. The time-dependent quantity (V) pinpoints the actual location of the 
satellite on the orbital path. These elements are (see fig. 3.1): 


a (semi-major axis) - defines size of the elliptical orbit. 
e (eccentricity) - defines shape (non circularity) of orbit. 





1 (anclination) - defines the angle between equatorial plane normal (K) and 
orbital plane normal (W). 


€? (longitude of the ascending node) - defines the angle between the vernal 


equinox direction (I) and the point where the satellite crosses 
the equatorial plane from southern to northern hemisphere 
(ascending node). 


( (argument of perigee) - defines the angle between the ascending node and 
the satellite's closest point of approach (perigee). 


v (true anomaly) - defines the angle between perigee and the satellite's actual 
position. 


The elliptical orbital path, is described by the equation of a conic section 


" a(l-e”) 


= (3.1) 
l+ecosv 


where r is the magnitude of the satellite's position vector at a specified v. In 


general, v is nonlinear with respect to tirae, so the mean anomaly M, which varies 


linearly with respect to time for purely Keplerian motion, will be used as a filter 


state instead. M is related to v through a quantity called eccentric anomaly E. It is 


defined as follows (see fig. 3.2): 


MzE-esinE (3.2) 
oa od (3.3) 
]-ecosE 


"S sin EN1- e? 


(3.4) 
]-ecosE 


Eqs. 3.1-3.4 define the relationship between M andr. 


The filter state vector will be 


CC A MJ (3.5) 


where 
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and Ty, is the time between observations. For the classical 2-body case, it can be 
seen from eq. 3.6 that state propagation is a linear problem. 

In order to improve the estimate of the state vector using available 
observations, 1t becomes necessary to compare the states with the observables in a 
common reference frame. Satellite position, which is implicitly a function of the 


states through the relationship in eqs. 3.1-3.4, can be written 


f COS 


Pis Gy) 


0 


where rP denotes a position vector in the orbital (PQW) reference frame (see fig. 
3.1). The remaining states are used to transform satellite position from orbital to 


inertial (IJK) coordinates using 


A 


Figure 3.2. Eccentric vs. true anomaly. 


cmcQ2—swsQci -sac -cosci sOQOsl 
% l , s 
C ^ -|cosQ-csocQci -sæsQ+caxļQci -cQsi (3.8) 
YANI CWSI el 


where c8=cos8 and sO=sinO for brevity. (See App. A for development of 


CA) 


Then satellite position in inertial coordinates will be 


T COS V 
p d rsin V (3.9) 
0 


Now. inertial satellite coordinates need to be related to a measurement in the site- 


based reference frame. 


B. OBSERVATIONS 

The detection scheme modelea here 1s a fixed planar radar beam generated by 
multiple continuous wave transmitters along a great circle path (see fig. 3.3). Asa 
satellite penetrates this beam, phase information 1s collected at one or more of the 
receivers. Antenna geometry at each site makes it possible to compute direction 
cosines, which are the pseudo-observations available for element set improvement. 


These observables are (see fig. 3.4) 


mee 
P, 
(3.10) 
cosß= EE 
Ip, 


where East-West cosine = cos a and North-South cosine » cos D . (y can be 


calculated from: cos? y 4 cos? B cos? & 2 1.) 





equato! 


Figure 3.3. Radar plane geometry. 


We have 
pcosa —p- E 
(3.11) 
pcosB=p-N 
and 
Ro pone: 
=| | (3.12) 
cosß 
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N (North) 





Figure 3.4. Receiver site geometry. 


is the vector of observables to be used for orbital parameter improvement. (See 
App. B for development of Z in terms of the orbital elements.) 


The position vector in the site-based frame 1s 


r" =CAr! (3.13) 
where 
cleo clso sf 
CA = —sodc(az)—sfcds(az)  c(az)eb —s(az)s£s X s(az)c4 (Sa 
s(az)sp— c(az)stcdb -s(az)cb-c(az)s{sb c(az)c/ 
and 


£ = latitude of receiver 


OE Elon. 


1] 


W = earth's rotation rate 
to = time of last Greenwich crossing of vernal equinox 


lon; ~ longitud o TEENE 
(Seehe >) 


(See App. A for development of E) 


Recapitulating, 


x, I ae) 


z -g(x,.T,) (3.15) 


both of which are, in general, nonlinear functions of the orbital elements x,. 
Assuming for now that perturbations to the two body dynamics and observation 
noise can be modeled as additive zero mean, white noise processes, the state and 
measurement equations become 

X A 


MENOR 


va 
(3.10) 
Z, = gx; TM 
Kalman filter equations are then 


A ep T dE ES Tea | 


re RT 
Pye = PP, ® +Q 


G =P, TP AT+R] 
E Ps 


P, -(I-G, H)P,,. (3.17) 






I EE. 


equinox 


Figure 3.5. Receiver site geometrv. 
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where 





à = ATi) 
Ox. 

é T 
ii = dg(x,.T, } 





(See App. B for derivations of ® and H .) 


(3.189 


In words, the actual observations, z,, are compared with computed 


observations, Ly (based on the classical two-body propagation model), resulting 


in the Kalman filter innovations. These are multiplied by the filter gain, which is 


based on the dynamic model and measurement uncertainties and used to produce 


an improved estimate of the orbital elements. 
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C. RADAR PLANE INTERSECTION 

Meaningful application of the observations to orbital element set improvement 
requires the ability to predict subsequent penetrations of the radar plane bv the 
satellite in question. A simple way to do this is to define the satellite's position 
vector in a coordinate frame referenced to the radar plane (See fig. 3.6). Then the 
out of plane component (Z) of position can be checked for a value of zero, 
yielding an in-plane condition. 

Satellite position, readily available in orbital frame coordinates (rP) is 
transformed to radar plane coordinates by performing two coordinate 
transformations: one from orbital to inertial frame, the other from inertial to radar 


plane frame, or 
DC (3.19) 


where r* denotes a position vector in the radar plane (A YZ) reference frame and 


co, co co, sê soe 
C% =| -s0 co 0 (3.20) 
S ROUES eye 


where 
EE (teats, | long 


0,,1s the inclination of the radar plane to the equatorial plane 
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lon, is the longitude where the maximum latitude of the radar plane occurs, 
which describes the location of the X-axis of the radar plane 
reference frame. (Scenes ch 


(See App. A for development of CA) 





Figure 3.6. Radar plane orientation. 


Since the only component of satellite position of interest in determining radar 


plane intersection is the Z-component, the expression can be further simplified to 
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AG 1]C4C%rº 
= [-s0 „P,c8 - s8 „P488 + c8 „P; |r cos v CoA. 
+|-s8,,P,c8 — sð Ps — c6, P, |r sin v 


where 
E > = 
a at (3.22) 
Es EG 


The in-plane condition of relevance to the filter is that which falls inside the field 
of view of one of multiple receiver sites. also placed along the great circle path 
circumscribed bv the transmitters. 'The time epoch of this condition is the time 
close to which we expect to observe the satellite from one or more of the radar 


plane receivers. 
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IV. DOMINANT PERIURBATIONS 


Satellite orbital parameters deviate from Keplerian motion deterministically 
and randomly due to many factors. Ignoring for now such intentional human- 
caused phenomenon as plane change and stationkeeping maneuvers, the largest 
natural contributor to parameter vanations 1s the nonuniformity of the gravitational 
field due to earth's oblateness. For low earth orbits, the next most significant 


factor is atmospheric drag. 


A. GRAVITY PERTURBATIONS 

The largest contributor to perturbations to the Keplerian orbital elements 
arises from the oblateness of the earth. These perturbations take the form of a 
combination of secular and periodic variations of the classical orbital elements. 
Many solutions have been developed for this problem, each with its own sirengths 
and weaknesses. All the solutions, however, are based on the gravitational 
potential 


u= Be 1-¥ 1 5,P,(sind)} (4.1) 


T n=2 D 


where 
N. ls earth's gravitational parameter 


r is the distance of the satellite from earth's center 
J, are n™ order zonal harmonics 
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P,,(sin ò) are Legendre polynomials 

Ò is satellite declination 
and are more or less complex depending on the order to which the terms of the 
expression are carried. This expression for U already carries with it a degree of 
simplification in that it assumes axial symmetry for the earth. In truth, the equator 
is slightly elliptical, but the main effect is on the stability of geosynchronous 
satellites, since their orbits are in the equatorial plane [3]. Further simplification is 
accomplished by dropping terms for n 2 3. This is justified for medium accuracy 


applications because of relative magnitudes of the J, terms 


IE Moe IG 


mc 2302 
(4.2) 
IN o 
where J, successively decrease. So the simplified potential function is 
u=(He+ Bes, (sin? 8-1)] (4.3) 
DOT 


Defining R= iG sin? 6-1) as the perturbing function, the variation of that 
I 


function with respect to each of the orbital elements is [ 4 ] 
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OR__3p 
oa a 


OR Y 
de 2(1-e?)r 


-sin i -sin 2u -sin v(2 + e cos v)] 


[cos œ(l+ ecos v)(1-3sin? i-sin? u) 








A dt 

A ASCO SA 

ol n 

R -3J, .,. . 

DN = sin’ i-sin 2u 

OQ 2r an 
R o 

O 

OR er dix 

OM n dt 


To develop a first order theory, the right-hand members of eq. 4.4 are 
expanded in a Tavlor series about the initial orbital element values, retaining only 
the first term of the expansion, yielding expressions fc; the derivative of each 


orbital element with respect to true anomaly v. Integrating these expressions 


vields an analytical first order solution in the form 


“Ed TO AVOA N) (4.5) 


where 
x = instantaneous orbital elements 


Xo = initial orbital elements 


Ó,x — secular variations as a function of v 


ôx = periodic variations as a function of v 
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Al! six orbital elements experience periodic variations. However, only Q2, @ 


and M experience secular perturbations. These are 


3,0 = E v.) 


p 





(4.6) 








5 (1-3sin* 6,)| (t-t,) 
r 


Note that 6, M is expressed in terms of independent variable t, which is the time at 
which true anomaly is equal to the value of vin òQ and ò œw. The periodic 
variations in €) and M do carry l/e, terms which approach singularity for near 
circular orbits. This can be explained by the fact that €) and M are undefined for 
circular orbits and ill-defined for near circular orbits. Fortunately, the combination 
of © + M obviates the singular condition by effectively forming one new orbital 


element out of o and M. 


Long term propagation of the orbital elements, assuming a drag-free 
environment, is accomplished by applying only the secular variations, since 
periodic variations will have no net effect on the elements over time. However, 
the periodic variations will be seen later to come into play when observations are 
assimilated into the orbital element improvement process. Inclusion of these 


secular variations in the dynamic model yields 


2] 


(4.7) 














x... =f(x,.T,) 
where 
a, 
Ct 
E 
E 3 
02, - 5 2 [cos(i, )(v.., =N ) 
3J, NS 
m te g 1 Cavs 
(4.8) 
3 % 
No E (8) (1-3sin 8). N 
2 
where 


p, =a, (1 —e; ) is the semi-latus rectum of the orbit 

J; is the second order zonal harmonic (equatorial bulge). 

The gravitational force model which gives rise to these secular variations is 
relatively simple. It assumes that the earth is symmetric about its spin axis and 


about the equatorial plane. Noteworthily, the first three mean orbital elements 


experience no secular variation due to earth's oblateness. 
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It may be possible to maintain identification of some satellites using secular 
variations alone. Continuous identification of a satellite, however, may depend 
upon the application of the periodic variations to the mean elements at observation 
times. These periodic variations are assimilated as an additive (superposition 
assumed) adjunct to the existing dynamic model and therefore are not technically a 
part of the propagation dynamics. The filter innovations resulting from a 
comparison of observables with predicted states are then applied, via the Kalman 
gains, to the mean elements, implying a prior removal of the periodic variations. 


These periodic variations,to first order, and ignoring terms of the order of e, and 

















smaller, are 
da, = STR 1. cos|2(o, + v) 
Be = sk cing. EE i,:cos(20, * v) 
MED 2 4 ° 
-— i, :cos(20, e$) 
12 
Jo 
ôi, =- " sin(2i,) cos|2(,  v)] 
3J 
00, = 7 —cosi, -sin[2(0, + v) 
3J 
00, = A = cosi, -sin|2(w, ar v) (4.9) 
3J 
ôM =- "m t -> sin? i, Jsin[2(o. +v)| 
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B. DRAG PERTURBATIONS 
Satellites traveling in the altitude regime below 600 km experience non- 
negligible drag force due to atmospheric resistance [3]. This force is well 


represented as 


] A 
i, = 7 00 (4.10) 


where 
m=mass of satellite 
Cp=2 1s a dimensionless drag coefficient 
A= effective cross sectional area of satellite 
V= velocity of satellite with respect to atmosphere 


p=local atmospheric density. 


1. Atmospheric Drag Effects 


To first order, drag force doesn't affect i and €) [6]. The main effects of 


this drag force are to bring about secular decreases in a and e, 1.e., to decrease the 


size and eccentricity of the orbit. @ and M also experience changes due to drag, 


but these are negligible compared to the oblateness effects. 


It turns out that 


(4.11) 


CLA 72 
da EI (1-0) 
dM m 


l-ecosE 


where 


di] —e costs) 
l+ecosE 


0 A) : 
dz " (l-e ) COSi (4.12) 





and 





Y 
de | Co^ i-e (1+ecos E)” 


- (1-1) cos E 
dM m (creas Ey? 


(4.13) 





E 
2m \l+ecosE 


7 ^ 
(1- t)sin E 


All quantities in these two expressions are well known (according to the simplified 
model) except for atmospheric density p, the value of which 1s related to space and 
time in a complex manner. Atmospheric models delivering as good as 15% 
standard deviation are in the form of tables or very complicated empirical 
formulae or both [13]. For simplicity, it is advantageous to use an analytical 
function to model p. It 1s assumed that variations in p in the upper atmosphere 
(150 km to 750 km) are a superposition of four individual variations [4]. 

The first variation is diurnal, or, daily. This fluctuation has to do with 
change in solar flux intensity associated with moving from darkness to daylight, 
and vice versa. Density, depending on latitude and altitude, may change by up toa 


factor of 10. 
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The second type is related to the 11 year solar cycle which is motivated 
by intensity of sunspot activitv. Maximum density, occurring at the peak of 
sunspot activity, is roughly ten times the mean density. 

A third variation is semiannual, with density falling to a minimum in July 
and a maximum in October, with a less pronounced minimum and maximum in 
January and April, respectively. The range of this variation is, generally, less than 
1/20 the range of the 11 year solar cycle. 

The fourth type of variation 1s nonperiodic and unpredictable and is tied 
to geomagnetic storms rising from solar flare eruptions. Atmospheric density may 
increase by a factor of 10 from the geomagnetic quiescent value. 

Over the short term (days), barring magnetic storms, the diurnal vanation 
dominates. However, over several years, the solar cycle has a more profound 
effect on atmospheric density. 

The modeling approach followed here was to start with a very simple, 
even crude, atmospheric density model, including the above considerations as 


necessary to achieve desired accuracy. 


2. Simplified Drag Effects 
First, simpler expressions for changes in a and e due to drag can be 
developed by viewing these changes as a result of orbital specific energy being 
extracted by drag force [14]. Specific energy of an orbit, which is conserved for 


an unperturbed Keplerian orbit, can be expressed as 


Bes “a 
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The change in a for an incremental change in € is then 


Da 


da= de (CDS 
" 





If it is assumed that most of the drag occurs near perigee (atmospheric density 


exponentially decreasing with altitude), then 


r=r, =a(l-e) (4.16) 
from eq. 3.1, and 
^ T 
iac m 
i B (4.17) 
- u(1* €) 


from eqs. 4.13 and 4.15. For a satellite moving through dv near perigee (see fig. 


4.1), the specific energy extracted will be 





Figure 4.1. Satellite moving through perigee. 
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F 
dc PM 


o (4.18) 
CA s 
2\ m 
Substituting eq. 4.17 into eq. 4.18 yields 
Li eae 
de=-- put +e)dv (4.19) 
2 


Substituting eq. 4.19 into eq. 4.15 gives an expression for decrease in semi-major 


axis due to a decrease in orbital energy 





em 
da = 2 P (1+e)pdv (4.20) 
m 
Since the new orbit will pass through the same perigee point, all of da shows up as 


a decrease In apogee radius, and 


dr, 7 da(1- e) - ade 2 0 (4.21) 


by differentiating eq. 4.16. Solving for de gives 


C E | 
der da (l-e)= 1 2 x za )pdv (4.22) 
m 


a 





3. Atmospheric Drag Characterization 


A crude model for variation of p starts with the assumption that 


temperature and chemical composition of a gas remain constant [14]. Then, 


density 1s directly proportional to pressure according to 
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p= Jr (4.2) 


(Realistically, k is a function of temperature and chemical composition.) At an 


altitude h, the decrease in pressure accompanying an increase in altitude is 
dP = —pedh (4.24) 


where g is acceleration of gravity. Since k is assumed constant, 


| 
dP 2 — dp (4.25) 
k 
Combining eqs 4.24 and 4.25 yields 
d 
cP = -kedh 
p 
P d h 
[ — - -kg| dh (4.26) 
o, P 0 
Pipe a 


So, In general, atmospheric density decreases exponentially with altitude. 


Eq. 4.26 can be rewritten 
ya p.e 7*. (4.27) 


where p, (reference air density) and h, (density scale height) change for different 


altitude regimes. For example, using standard fixed values of p from the "U.S. 
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Standard Atmosphere, 1976," a reasonable fit can be made by applving the values 


in Table 4.1 for the ranges of altitude shown. 


Table 4.1. DRAG PARAMETERS. 


|. BOkmezel0km  | ——— 185 | 58 
20kms | 12099 — | Ho 
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V. THE RECURSIVE ORBITAL PARAMETER ESTIMATOR 


An extended Kalman filter, designed for orbital parameter estimation, is based 
on the earth-fixed radar plane detection scenario described in Chap. III (Also see 
fig. 5.1). The filter propagates the orbital elements in time according to a first 
order nonlinear dynamic model which takes into account the dominant 
perturbations arising from earth oblateness and atmospheric drag. Upon receiving 
observations, in the form of pairs of direction cosines, the filter calculates an 
improvement to the current parameter estimate. The filter operates recursively as 
observations become available. The filter is written in Matlab code and can be 


found in App. C. 


A. FILTER INITIALIZATION 

The estimator is designed to receive two data files. One contains an initial set 
of orbital elements along with subsequent pairs of available direction cosines (see 
App. C). The other holds information about the receiver sites. Values of pertinent 
constants are also set, as are the initial times necessary for orbit propagation and 


estimation of observables. 
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INPUT DATA 


INITIALIZE KALMAN FILTER 
PROPAGATE STATES 






ESTIMATE MEASUREMENTS 
CALCULATE SP OE 
APPLY STATE IMPROVEMENTS 


NO 


Figure 5.1. Estimator flow chart. 


1. Input Data Files 
Two input files are required for filter operation. One contains an initial 
set of orbital parameters, along with subsequent sets of direction cosines. This 
information is all tagged with corresponding times (see App. C). The last five 
orbital elements occur in classical orbital element form, but, in lieu of semi-major 


axis, orbital period is given. The conversion is simply 


y 
Tu? 


DEE 
T 5 
(E 


where Ta: is the period of the satellite. Note here the apparent disappearance of 


earth's gravitational parameter ue. This convenience is brought about by the use 


of canonical units [1], where 
DU 
=] = 
me TUE 


1 DU¿=6378.135 km and 1s earth's mean equatorial radius. 





1 TU,= 13.44686 min and is the time it takes for a satellite traveling in a 
circular orbit at a radius of 1 DU to travel 1 DU of arc length. 


Canonical units are used in all filter calculations. The other input file contains 
information on the receiver site locations,altitudes, and azimuths. The azimuthal 
measurement indicates, for each receiver, the angular offset of its coordinate frame 


from true north. 


ES 


2. Constants 


The earth's apparent rotation rate is initialized as 


We = 05883378 ads, (5.2) 


The radar plane's inclination to the equatorial plane and the longitude of its center 


line are, respectively, 


O not: 
i (55) 
lon, = —101.31348° 
The radar plane is oftset from earth's center by about 20 km, or 
Loo e SH DE (5.4) 


3. Time 


All times are given in Universal Coo: iinated Time (UCT) and converted 


to TUes for filter use. Temwch. Which is time elapsed since the Greenwich meridian 
last passed the geocentric inertial I axis, is computed: from a reference angle 


(measured from I to Greenwich) at O hrs, 1 Jan 92, of 


0 = 98.920250931' (5.2) 


gwch 


using the @, given previously. 
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4. Filter Matrices 


The error covariance matrix P is initialized to reflect the level of 
confidence placed in the accuracy of the initial orbital element set. The 


measurement error covariance matrix R is set to 
aie. ai 
= > o) 
0 


where O... = Ic, -.0002' reflects the uncertainty of the available direction 


cosines. 


B. STATE PROPAGATION 
The filter states are propagated according to the two-body Keplerian orbital 
dynamics perturbed by earth oblateness and atmospheric drag. The nonlinear state 


equations are 





C- -AN 
d, | - jose. =v,) 
In 


e, (El ha -e’\p[v,.ı-V,) 





m 
l, 


[cosíi, Jv... E 


- 2-3 sin’ à f. -vi) (5.7) 


EAE e 
ves rela i-e) T, 


I, 


3J 
de e Q, -— 








ES 


T N 
and are propagated at arbitrary increments of. "7e . This is continued until a 


crossing of the radar plane is detected by a sign change of the radar-plane-normal 


component of satellite position. 


C. RADAR PLANE INTERSECTION 


Detection of a radar plane crossing calls a function which iterates to an in- 


plane condition for estimated satellite position. This 1s accomplished using a 


Newton-Rapson iteration on the expression for the plane-normal component of 


satellite position. The iteration calls for a calculation of an approximate AM 


which will drive rz to near zero. This is given by 


where mom ear dir 


and 


ANE — 


T, 
uz 


I = Arcos VEA TSn 





or, OX, a o(r cos v) 
OM OM OM 
o in V 
+ iin vic 
oM 


Taking the individual derivatives yields 
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(5.6) 


OF 


(5.10) 


where 


and 


where 











OX 0 
— — —- 
OM OM 
o 0 
: mg E ) ps 0 
OM "o 0M 
o(cos 0) N 
= -w.a 72(Sin 8 
EM oa 2 (sin 6) 
o(sin 0) -% 
= (0.2 (cos8 
N EN ova) 
cose — asinE 
OM l-ecosE 
Om ia EE GE 
OM ]-ecosE 


OE | 





fees = a\cosk—c) 


rsinV=avl—e’ sinE 


] 


oM — l-ecosE 





SH 


(5.13) 


(5.14) 


This iteration is found to converge to within z, <€ in three to four iterations, 


where e = 107°. 
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The function then determines whether the satellite 1s within the field of view 
of one or more receiver sites by calculating the in-plane r, and r, components of 
satellite position and checking the inequality (see fig. 5.2) 


r 

Sr (52159 
I 

It should be noted here that the field of view needed to be opened up to pick up 


some satellite observations by receivers at the edges of the radar plane. 


within field À 


MEN. ot 
l (lon « 101^) 





Figure 5.2. Sensor field of view. 


D. ASSIMILATION OF OBSERVATIONS 
Once a Satellite is in the field of view of radar plane coverage, the input data 


file is checked for observations at or near its time of arrival. If a pair of direction 
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cosines is available, then a measurement estimate is calculated according to the 


> 

Z = Su 

Z A (5:16) 
p 


® is also calculated according to App. B. The plant noise covariance matrix Q is 


development in App. B, where 


calculated using the sguares of expected levels of uncertainty in each of the states. 


For example, if the two-body problem is the dynamic model, then the Q matrix 


would be 
o: 50 0 0 0 0 
Ob tuo 0 0 0 0 
0 0 0° 0 0 0 
Q= a (ED 
oe oes “6 0 
0 0 0 pos “O 
0 0 0 0 D "Gn, 


where c, are proportional to the magnitude of the largest unmodeled perturbation, 


in the two-body case, earth oblateness. So, for the two-body problem, 


py 


Eua 





(1-e*)sin* i-cos[2(@ + v) 











2p' 
35 
EE eos ve zsin! iccoso e v] 
2p’ 2 4 
2. 3) 
o —= on Dr ee 2 en 
O, 2 8 2 | ( ) 
Gg T 
CM > 
era eu. (3.19 
2 n 2-3 sin ir, 
a 2 





3 Y 
LUE Je -e?)^ T, | 
E 


Now, the error covariance matrix can be calculated as 


P,_ =®P.., o «Q (5.19) 


/k-] 


H is then calculated as given in App. B, making possible the calculation of filter 


gain 


HT R) | (5.20) 


which is then applied to the residuals z-2Z to calculate improved orbital 


parameters Xx. 


It should be noted that, throughout the filter, the quantities v and E are needed 


for various calculations. E is calculated fom M using Newton-Rapson iteration on 
M=E-esınE (5.219 


and v is calculated by solving for 
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j — 2 1 
V == sın =) v1-e' sinE x 3,22) 
1-ecosE 
then using 
E-e 
ve cos EI (9-95) 
l-ecosE 


to resolve the ambiguity. 
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VI. RESULTS 


This extended Kalman filter is designed to estimate the orbital parameters of a 
satellite, given observables in the form of a pair of direction cosines obtained upon 
the satellite's passage of an earth-fixed, stationary planar radar beam. One day's 
worth of data for each of three satellites was obtained for filter testing and 
adjustment. Results of the application of three successively more comprehensive 
filter models are compared with each other. 

The first filter models only two-body motion, the second adds secular 
perturbations due to earth oblateness, and the third also includes periodic 
perturbations due to oblateness as well as atmospheric drag effects. Each filter, of 
course, models unknown perturbations as zero mean, white, Gaussian noise, in the 
tradition of Kalman filtering. Unfortunately, none of the available satellite data 
include an orbit in the drag regime. The filters will be compared in the areas of 
radar plane time of arrival, a residual small sample statistic, and orbital parameter 


estimation. 


A. RADAR PLANE TIME OF ARRIVAL 

With orbiting objects arriving at the radar plane on the order of every few 
seconds, it is important that prediction of arrival times be fairly accurate. Current 
methods yield accuracies of within a second. The best that can be hoped for under 


the current detection configuration is about 0.25 second, which is the sensor time 


uncertainty. After some crude tuning of the three filters, a comparison of arrival 
times for the three satellites is shown in Table 6.1, where At - t, —t... 

It can be seen that, in all cases and for all models, the first detection time after 
the estimator 1s started up (which occurs on the order of 12 hours after filter 
initiation) is always mispredicted by the filter by a time on the order of seconds. 
However, it is seen that the secular and periodic models cut down that time error 
by a factor of from three to four. Subsequent detections yield differences in time 
between calculated and observed plane intersections of less than 1 second for a 
one-orbit-later detection and between two and three seconds if detection occurs 


next on the "backside" of the orbit (about a half day later). Although this is true 


Table 6.1. OBSERVED VS. ESTIMATED TIMES OF ARRIVAL. 


Satellite 116 


Satellite 117 


Satellite 118 
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for the two-body and secular cases, the incorporation of periodic perturbations 


yields subsequent time differences of less than one second in all cases. 


B. RESIDUAL SMALL SAMPLE STATISTIC 
The filter residual is the difference between actual and estimated observations, 
and, as such, is the basis for a good performance measure. The small sample 


variance of the residual 1s calculated for each case as follows: 


necs Lm 
c 2$ —-X' (6.1) 
1z1 n 
where 
a Ex 
xX=)- (6.2) 
iste 


Is the small sample mean, x; being the 1" residual. Plots of these variances for each 
case are shown In figs. 6.1 - 6.3. 

It is evident that the residuals for the N-S direction cosines are always very 
small. This is because this direction cosine is always close to 90° and the 
geometry of the problem always forces the estimate of this direction cosine to be 
very close to 90°. However, the E-W (in-plane) cosines are more interesting to 
observe. In general, the trend is toward smaller residual variances from 
observation to observation and from simple to more complex dynamic model. The 
only case where this does not hold true is for Satellite #117, where it can be seen 
that the two-body model is sufficiently lost after approximately eight hours to 


yield higher residuals than at the previous set of observations. Note that a "stack" 


of asterisks denotes several observations gathered at the same radar plane crossing. 
which are then applied one at a time for parameter improvement. 
It can be concluded from this limited residual data that, for the relatively short 


times included in the data, the filter is performing acceptably. 


C. ORBITAL PARAMETER ESTIMATION 

The ultimate objective of this filter is to be able to predict what a satellite's 
position will be at some future time. Hopefully, the filter 1s predicting plane 
intersection times bv accurately estimating orbital parameters. Figs. 6.4 - 6.6 
show time propagating values of the four parameters that are expected to be 
varying most with respect to their two-body propagation values. All three 
satellites show similar results. In all cases, semi-major axis estimates differ 
greatly from model to model. Eccentncity estimates for the two-body and secular 
models are very similar, with the periodic model showing the greatest deviation. 
All estimates are very similar for the longitude of the ascending node, though, 
between observations, the two-body model does not have a mechanism to 
accurately propagate this element. Argument of perigee estimates are very close 
between the secular and periodic models, while the two-body makes very little 


correction to this parameter. 
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Figure 6.1. Variance of residuals (sat. # 116). 
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Figure 6.2. Variance of residuals (sat. # 117). 
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Figure 6.3. Variance of residuals (sat. # 118). 
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Figure 6.4. Orbital element estimates (sat. # 116). 
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Figure 6.5. Orbital element estimates (sat. # 117). 
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Figure 6.6. Orbital element estimates (sat. # 118). 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


Orbital parameter estimation was attempted using an extended Kalman filter 
for each of three orbital dynamic models. Observations were angles only and 
made roughly twice a day per satellite passing through the radar plane. It was 
found that the two-body orbital model, including secular and periodic variations 
due to earth oblateness, was able to predict subsequent radar plane crossings to 
within a second of actual crossings. 

The filter, which includes a crude model for atmospheric drag effects, was not 
tested on anv satellites in the drag regime. Future research should include data 
from satellites with altitudes less than 600 km. Also, longer term data should be 
obtained to check filter performance over the long haul and also to see if the filter 
could track longer term disturbances or perturbations to the orbit under 
consideration. 

Tuning of the filter could also be accomplished by obtaining a larger data base 
of satellite orbital elements and corresponding observations. This could be done 
manually, but it would be more exciting to See if an artificial neural network could 


be designed and trained to accomplish the same task. 
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APPENDIX A 


COORDINATE TRANSFORMATIONS 


The calculations and algorithms presented 1n this paper necessitate describing 


vectors in several different coordinate frames. These frames are all orthogonal. 


Transformation of a vector from one reference frame's coordinates to another's is 


accomplished using direction cosine matrices (DCMs). 


Reference Frames 


The reference frames of interest are as follows (see fig. A.1): 


Geocentric inertial (1JK) - The I vector points inertially in the vernal equinox 
direction, with I and J lying in the equatorial plane, and K 
piercing the north pole (coincident with earth's spin axis). 

Orbit-fixed (PQW) - The orbital plane is the fundamental plane, with P 
pointing to perigee (point of closest approach), Q 1s the semi- 
latus rectum direction, and W is the orbit-normal (also orbital 
angular momentum direction). 


Radar plane (XYZ) - The radar plane is the fundamental plane, with X 
pointing to the maximum north latitude point, Y lies in the 
equatorial plane, and Z is the radar plane normal. 


Site-based (HEN) - H 1s local vertical, E 1s east, and N is north. The origin 
coincides with the radar site location of interest. 


Direction Cosine Matrices 


Three DCMs are developed here, performing one orthogonal rotation at a 


time, then combining the rotations. 
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LIK (Geocentmc inertial) PQW (Orbit-fixed) 
K 


XYZ (Radar plane) HEN (Site-based) 





Figure A.1. Reference frames. 


Orbit-fixed to Inertial 
Satellite position 1s most simply expressed in orbit-fixed coordinates, 


but for comparison with earth-fixed quantities, it is useful to transform both to an 
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intermediate reference frame. The inertial coordinate frame (in this case, 
geocentric inertial) is the obvious choice. Fig. A.2 shows the individual rotations 


required to transform a vector from orbit-fixed to inertial coordinates. 


COS) Sin Q2 
rotl=|-sın cos 
0 0 


=> 


0 0 
COSI sin] 


—SIN1 cos] 


COS (D sin W 


rot22|-sino  cosQ 
0 0 





Figure A.2. Orbit-fixed to inertial rotations. 
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Combining the rotations yields 


C - [rot3][rot2][rot1] 


cmc Q - sosOci cwsQ2+smcOQci  SOSI 


=| -swcO —- cwsS9Qci —swsí2+ctcóc cosi N 
sQsi —sicQ2 cl 
Because the transformation is orthogonal, 
cwcS — swsSci -secGQ -casci sOsi 
c% =[C%] = cmsQ+sacQci —-s@sQ+cemcQei —sicQ 
swsi COSI ci 
(A.2) 
Bo Ry ak. 
EE OE EP 
Pep 


Inertial to Site-based 
Fig. A.3 shows the individual rotations required to transform a 
vectors coordinates from the inertial to a site-based reference frame. Rotation #1 
keeps track of earth rotation. Rotation #3 points out the fact that the E and N 
vectors are not truly east- and north-pointing. E is in fact tangent to the great 
circle circumscribed by the radar plane, and N is orthogonal to E in the local 


horizontal plane. 
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cos sinó O 
H rotis|—sind cosó O 
Ò 0 0 ] 
] 
where 0= Oe(t-t,)- lon, 
O; is earth's rotation rate 


t, Is time of last Greenwich crossing of I axis 


lon,, 1s longitude of radar site 


K 
cost O sin / 
H TOL = 0 ] 0 
f -sın{ O cos/ | 
I 
where Y 1s latitude of radar site 
K 
1 0 0 
E rot32|0  cos(az)  sin(az) 


Lom O —sin(az) cos(az) 
J 


Figure A.3. Inertial to site-based rotations. 


SE 


Combining the rotations yields 
C^ z [rot3][rot2 [roti] 


(A.3) 
c£có c£só Su 


= | -sobc(az)-sichs(az)  cóc(az)-s(az)s/sÓ X s(az)c4 
sÓs(az)— s£cóc(az)  —cós(az)— s£sÓc(az) c(az)c4 


Inertial to Radar Plane 
The transformation from inertial to radar plane coordinates is very 
useful to the determination of satellite crossings of the radar plane. Rotation £1 


accounts for earth rotation. Fig. A.4 shows the individual rotations necessary to 


perform this transformation. 


Combining these rotations vields 


| cd „cd cd „sd scm 
C% =| —sê cO 0 (A.4) 
SUED SUME) Com 
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cos sinO O 


rotl=|-sinð coso © 
0 O l 


where 6 = We(t —t,)+ lon, 
a% is earth's rotation rate 
t, is time of last Greenwich crossing of I axis 


lon, 1s West longitude of maximum North 
latitude of radar plane 


cos, 0 sinO 
X TO = 0 ] 0 


Er SAO O cost 
I 


where 6, is inclination of radar plane to equatorial plane 





Figure A.4. Inertial to radar plane rotations. 
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APPENDIX B 


LINEARIZED FILTER MATRICES d AND H 

The extended Kalman filter can be viewed as a predictor-corrector type 
estimator. States are predicted using non-linear plant dynamics. Then a correction 
is computed by applying the Kalman gain to the filter residuals (observations 
minus estimates). Calculation of the Kalman gain requires some type of 
linearization of the plant and measurement dynamics (see eq. 3.18). Normally a 
Taylor series expansion is used, keeping only first order terms which are evaluated 
at the best current state estimate. The general form of the linearization of a 


nonlinear discrete function 





Cx 
[x n 
si IE 8 .) (B.1) 
FE (x. By 
1S 
of, of, o 
Ox, OX, OX, 
of, : . 
(DEO S l (BZ) 
of. as 
Ox, OX n, 
brum 
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Development of Linear Plant Matrix ® 
As will be seen, the state matrix can be as complex as one lets it 


become. Following is a presentation of three progressively more complex models; 
1) Keplerian two-body motion, 
2) inclusion of dominant first order earth oblateness gravity effects, and 
3) dominant atmospheric drag effects. 


Two-Body Problem 


State dynamics for a Keplerian 2-body orbit are given as 


dis d, 
Cu. e, 
E à = f(x,.T,) (B.3) 
Qui Qu 
| Di. | OD 


Linearizing gives 


| ed io 
0 10000 

! 0 Oi O 

qe 
0 ONE oo dam 
0 OOO] CO 


3E 
2 k-1 


i 
| 
£5 
f 
Se 
E 
> 
= 
= 
em 
bad 
le a 
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Secular Variations Due to Earth's Oblateness 
To first order, the last three orbital elements are affected by secular 


variations, so that 











3 


n 


A 


E ma no 3 
M após “(1-35nº6,)] T, 


— 


where 


= Pi 


ja SES 
Vie. oro 


(B.6) 


sind, =sini, sin(v, +0,) 
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Linearizing gives 


























] 0 0 0 0 O 
0 ] 0 0 0 0 
0 O l 0 0 0 
s 
OD DO Qo ] 0 OO, 
ga, de, Qi, oM, 
om... oo ME ele o ' OO... (B.7) 
Gan de, di, OM, 
ONA or TEON o OBS NUNT 
Qa, de, di, AO, one 
where the terms in the fourth row are found to be 
eV) op. 
i =—cosi, AV P. 
Qa, p, da, 
Coa S| Op. 
= — cosi, AV dy 
de, p, e, 
AQ 3J 
— = sini, Av 
di, 2p, (B.8) 
OO. _3), m Qv, 








oM, 2p! < oM, 


65 


where 





pe 
da, 
9p, = —2ae 
de, (B.9) 


ov, I(I+e?)sinE+2esinEcosE 


OM -sin v, (l -e cos E)” 


k 


The fifth row terms are found to be 


where 


OO... -3J Q 
Es EN AS 
Qa, p,  * da, 











OQ) —3J, o 
K+] = er A, Av Pi 
9e, [ox NEM S 
dm, —I15J, .. . 
EE SITIS COST TAN | 
gi, 2p- (B.10) 


aa ale A OV, 


AM, 2p? “OM, 








A, = 2-2 sin? i (B 1 
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Sixth row terms can be evaluated by 


where 





























oM,., » — A, A, T, Or, 
Qa, D N TN 
M -9J or 
o k+l = - A, A. m k 
de, I NEUE 
OM -9J o(sin Ô 
un = — A, sin 6, T, ( * 
Qi, n : Qi, 
ON ON o(sin 6, | 
se A, sinó, T, ————— B.12) 
JO, P. 64 E do, ( 
oM, QA, Qv 
4 k +1 = T A, l 8 \ k 
oM, 2E en ON. 
As, 
— A. 


ale prcosv 


AT 
STIS 2. 00 O(sin 8 
ur een, Or Sin, 8. 
ord OV, OV, 


Or, 
OV, 





d(sin 6, ) 


OV, 


2 a 
a, (l-e; jel Sin V, 


WM CEN (B.13) 


(1+ecos v, )' 


- sini, cos(v, +0, ) 
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and 


A, -1-3sin'Ó, 


k 


J Y 
2 
2 (SEA) 


(B.14) 
A, =1+€, cosv, 





A, = - (1-3sin’ 8,) 
I, 


Development of z and H 
The relationship between the measurements (pseudo-observables) and 


orbital parameters (filter states) is contained in the nonlinear function 


O (B.15) 


and is the same regardless of the state dynamics model being considered. The 


observables of interest here are the north-south and east-west direction cosines, or 
cosa 
LE (B.16) 
cosD 


and are related to the range vector p by 


Pe 
a Ban 


68 


Note (from fig. B.1) that the HEN coordinate system in which. p is received is not 
truly up, east and north. H is local vertical, but E is tangent to the radar plane 
great circle, and N is orthogonal to E in the local horizontal plane. This site- 
dependent information is included in the transformation matrix CA (see App. A) 


as the anele az between true east and E. 





27; 


S QS 


Figure B.1. Site-based coordinate system. 


SO P;, Px, and range magnitude p need to be expressed in terms of the 


orbital parameters in order to form z = g(x,T). It can be stated that 


Pu p cos 
pg peosa (B.18) 
Px | LpcosB 
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In order to describe the direction cosines in terms of orbital parameters, first 
consider fig. B.2. Earth's flattening has been exaggerated in the figure to clearly 


show the effect on radar plane orientation. From the figure, it is clear that 
r=r,.. tR.. +p (BOSE 


where 


R.. - R(lat) * alt, 


site 


(B.20) 








Figure B.2. Relationship between satellite position and range. 
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R(lat) 1s the latitude dependent radius of a point on the reference geoid (due to 
flattening) [7] and is assumed to be in the H direction. Tofiser 1S assumed to be 
totally in the negative N direction, and alt... is the receiver's altitude. Combining 
the elements of eqs. B.18 and B.19 yields 

jc Deos 
r = pcosa (B.21) 


—T ay FPCosB 


rHis the transformation of r? from the orbit-fixed coordinate system to site-based 


coordinates, or 


p ey" 


ba ED ECOS 
O EP, TS ELS B 
LS ELI OE UD. P, 0 (B.22) 


P, 
E 
P, 

r(R, cos v+R, sin v) 

- | r(R, cos v R, sin v) 
r(R, cos v R, sin v) 


where 
es 
R,=H,P,+H,P,+H,P, 
R, =H,P, +H,P, +H,P, 
BEE EL HE ER P, (B.23) 
R, -H,P,-H,P, 4 H,P, 
R, =H,P, +H,P, +H,P, 
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Combining eqs. B.18,B.21, and B.22 yields 


pe r(R, cosv+R, sinv)-R 


PERS 


Pr 


site 
r(R, cosv+R, sinv) (B.24) 
(R, cos v R, sin v)- r,. 


From eq. B.24, p can be calculated as 


p=/P, +P: +Px (B.25) 


Now, eq. B.15 can be used as the estimate of the observation to compute filter 


residuals. But the Kalman gain depends on the matrix H, which 1s derived by 


linearizing g(x,, T, ) with respect to x;. 


The form of H is as follows 


oM (B.26) 
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where 


p 
(Po) ra 1d 


da p da =p? da 


p 
a Ni ] op 








M oM >M 


Dx 
a a OD. 1 op 








da p da s oy Gk 


Px 
a A) 100 ] op 


OE A Mi 








Following are the details of this development. 


op 
Oa 





Px 
Qa 


dp. 





9p 
Qa 


- (R, cos v & R, sin ao 
da 


=(R,cosv+R, sin yo 
oa 


- (R, cos v € R, sin JE 








+P; FPn 


da da 


p 
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+ OP 9p, 9p. 
2 Qa 


| 


(B.27) 


(B.28) 


where 


or la 


— = — B.29 
oa 1+ecosv ( 


The expressions for the change in p with respect to eccentricity e are identical to 


eq. B.28 except for the replacement of a with u which is 


or -a(2e+cosv+e” cosv 
de (I+ecosv) 


Since the R, terms in eq B.24 are functions of the orbital elements 1, (2, and @, 
the expressions for the change in p with respect to these three elements are as 
follows, with the replacement of ol by the appropriate element: 


a = = cos V+ e sin v} 











m - 


- - 
w - 








9p, OR, oR . 
=| —— cos V+ sinv T 
OC o< Oz 








OD. 
PN _ aie EES ed vlr 
OG Oc OC 


(B.31) 








op al Pu. OPs | 
H 


—— + jes 
da p ac gs Oc n dc 
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The expressions for the change in p with respect to a change in mean anomaly 


M is more involved because of the relationship between mean and true anomalies. 


This relationship is repeated here for convenience: 


























M=E-esinE (B.32) 
COS — ER de (B.33) 
¡Secos 
: / ced 
sin V = sinEyl=e" (B.34) 
l- ecos E 
Taking partial derivatives as appropriate yields 
Q O o(sin v 
erus ke ECT n DR r+[R, cos VR, sin yg 
oM OM OM OM 
o o In V 
P: IR AAN Pn r+{R, COSV t he sin js 
OM OM ^. oM É OM 
O E O j 7 , 
DBIS TR. (cos V) | R. HR. EOs VR, si v|— 
OM OM OM OM (B.35) 
m OD; ; ODE d OP. 
oM pl" aM “MM "aM 


18. 


where 


o(cos v) N ov OE 


dM OF OM 


A(Sin v) OV OE 


= cos V—— — 
dM JE ON] 





or ae(1 — e? )sin V ov OE 


oM (14e cos vy” dE ƏM 


ax -(1- e? )sin E 42e sin EcosE 


dE -sinv(l-ecos E) 
(B.36) 


aE 
OM  l-ecosE 


All the expressions for the partial derivatives of the Rn terms can be found in 


Tabie B-T. 
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Table B.1. DEVELOPMENT OF TERMS FOR H. 


OR, 
no Ossos + se, 


aR. 
E =P,H, "bL PH, 
Oo : 


- Ss si — s06O0s8H SE SOxTH. 
1 
OR, 
=-P,H, PENT 
oQ 


OR, 
=P.H,+P,H, +P,H, 
AO À 


R, 
E RAS GDS AS Tie SO DE SEL]. 


l 


OR, 
— = EH. Hobbs 
o | 


OR, 
m rer He op, 


OR 


1 


== cwsNsiH, - cweQsiH, + cwciH, 


OR. 
p — E H, +EH. 


oR 
— =-PH, -P,H, -P,H, 
Qo 


OR, 
Ee ed Did ees, Hoof, 
l 


OR, 
= Pil PH; 


OR, 
SP EH, zt 
dW 


— = cosQsiH, ^ coxcOsiH , + cwciH, 


1 


OR, 
= EE GP OH; 


OR 


—: --PH,-P,H,P,H, 


O 





APPENDIX C 


FILTER CODE (MATLAB) 

Following is the Matlab code for the extended Kalman filter described in this 
paper. This particular filter models the 2-body orbital dynamics problem, 
including secular and periodic perturbations. It requires 2 data files for input (see 


header of filter code): 


datain - contains a time-tagged initial orbital element set and an unspecified 
number of time-tagged pairs of direction cosines and 
corresponding receiving sites 


siteinfo - contains pertinent information about the set of receiver sites 
comprising the radar plane detection system. 


The filter calls 3 functions: 


enewton - performs a Newton-Rapson iteration to obtain eccentric anomaly, 
given a mean anomaly and eccentricity 


intper - performs a Newton-Rapson iteration to propagate satellite position to 
an ın-radar-plane condition. 


twopi - forces angular measurements to a number between 0 and 27. 


Extended Kalman Estimator 


The code for these functions follows the filter code. 
O0 * Sk Sk ak ak ak ak ak ak ak ak ak ak ak ad ak ak ak ak ak ak ak ak ak ak ak ak ak ək ak ak ak ak de de 3k N de ok oke ok ok de de 3k 3k a a ed ld al dl a 2k kc ke K 


This estimator requires 2 input files om 
DATAIN & SETEINFO, 
datain=[250 7116 82 19162602977 Il 
103.591322 .00756997 66.981264 1850 NNE 26.1 J 689 
12625310] 
83 65024.876 6 14583 972, 00750 
83 BSA36.2253 STE EA ENG 
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"63055 MM 
262952) eo 2 
-5 74050 E s 
SUED OD NE OE 
52. AOS O 
RR SA ROLO 
soe 405 107 
Ez .ol6u 0 1 
29290 RE O 


es 6343652529 
83 6455... 3 97 
83 63436.406 
83 03223500025 
o3 iese T2z 
e visada sa) 
83 174304.904 
eon 2743085900 
637174304906 


ON I OO AI od 
O O O OO COI 


So alza 340  =—11le. 97016 1.87473e-5 8.58120 
EA AO UOU24 4Z. 21244 9e-4 3.13900 
DES O CS COB DE q)  sA 21563 
Ee ops cT 56 9950I346—-7———5.66895 
EE NS. EE: AA PES -9,71924 
PA. 5] 9009 3 56002e06—565 10.5800] 


OY Cn iS CO h2 


MEM EA OE DR ARA EE EES EE RE EK KOEL ER E di 


ucueon[xkeepb,Tkeep, Ts mt,delTint,resid,xint,Pkeep]-... 
open idataln, siteinto):; 
m Se)ls-[datain(l,:)';datain(2,:)']; 
Be ry-dacaın (3: lenctchldetain(:,1)),1:5); 
K=inpels (1); 
ser 1=1:7 


Rar (a)=siteinio WEE SO: 


(zZ 
MI site ol), SS) *pI/180: 
mle (y= Sitei mao, 4) ; 
Be) sıreimto(ı1,5)*r91/180; 
end 


ARARAS AA RÃ A A 


DUC PROGRAM PROPAGATES A SATELLITE ORBIT AS A 2-BODY PROBLEM 
WITH SECULAR VARIATIONS IN OMEGA, omega, and MEAN ANOMALY. E 
Sees THE FOLLOWING FUNCTIONS: 

mYVOPT — places angle between 0 & 2*p1. 


Hi TEE - iterates to an in-fence-plane condition when a fence 
plane crossing is detected. 
ENEWTON - solves for eccentric anomaly, given values for mean 


anomaly and eccentricity. 


SEC kk ke ke ke RR AR RR RR RR RR RR RR RR RR RR RR RR RR RR RR e AR AR RR RA RR RR ER RR RR RR RR RR KOR RR 


E DEFINE EARTH ROTATION RATE * 
SS Ck C Ok Ck Ck Ck Ck C CC C Ce OI AAA AA A A AAA AA AA A A AAA AA A A A 
eearchrot=.05883378171654: % Define earth rotation rate (rad/TU). 
02=1.082645e-3; 


SARAH AAA AAA AAA RARA 4 


ME DEFINE VARIOUS TERMS + 


SS kk Ck ck Ok Ck Ok OO OO OK CO OK AAA ARA AR 


R-[4e-8 0;0 4e-8]; % Measurement noise matrix. 
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Pkkml=eye (6) *le-8; % Initialize error covarianco ncia 
Pkeep(:,1:6)-PkRene 

resao(:,1)- [0205 

E zeros (6,2); & Initialize Kalman gain. 
fince323.580310*932 arco & Define fence inclination. 
COSTINCE=COS (ELLAS 

Sın ine sin (uno 

lonx=-101.31348*p1/180; % Define location of Fênec =a 


GK KKK KR RR KK RR KK RK RR AAA IA A AAA AA AAA A AAA AAA A AH 


= INITIALIZE STATES (ORBITAL EME) į 
SK we Rk RR RR ARA ok a kk KK RR RO OK 
#now=[ (inpels (7) / (2*pi* 13944269317 NN ee 
inpels (8) ;inpels (9.12) EA cae, 
xkeep(:,1)=xnow; % Save original states. 
SRK KKK KKK KK KKK KKK RRR EK KKK KR KKK KKK A AA CA AAA AAA KH U KR KKK 
= CALCULATE. INITIAL ECCENTRICTAN ON Di x 
EST NK CK KOC OC OON OXON KON E Oe COR CO CIC Ok Ox OX On REKE E 


EA (1) =Enew:z0n (xnow(0), ¿Now (6) ES d a 


Sox X AA AAA AAA AAA AAA RK KR KKK KKK KK KK KK 


x FIND INITIAL É ane Fa 
SRK KK RR RRR ERK KKK KR KKK KERR KKK KKK KER KKK KK KKK EKER RR KKK KK RK KK KKK KK KK 
COSTA=(cos (BA(1))-xnow(2)) / (1-xnow (Z) “ees (Bae 
sinTA=sart (1-xnow(2)°2)*sin(EA(1))/ (1-xnow (2) *Ccos (2a IS 
rF=xn0w(1)*(1=xnow(2)*2)/11>2n0w (2) "cose 
TA (LISLWOPAIAI(ZSID (SINTA se 

qu Be gue 

lf cosTA«O 


DA 1-519 (079 * Perform quadrant Chee. 
end % for true anomaly. 
else 


1f cosTA<0 

RE EA 
end 
end 


GRR KKK KKK KKK KR KKK RK KR KK EK KR KR KKK KKK KKK KKK AKER KKK KKK KAKA KK KKK 


* CALCULATE TRANSFORMATION MATRIX (ORBITAL TO INERTIAL) A 


CARR RARA RARA AAA RR AAA AAA RÃ A 


cosperi=cos (xnow(5)); 


Sinperi=sin (snows); 
cosascnd=cos (xnow(4)); % Calculate Sines & cosines of 
Sinascnd=sin(xnow (4)); $ angular orbital elements. 


) 
) 
COSTLAESEOS (USAS) 

sinisme (AN 
P1scosperi*cosascnd-sinperi*sinasend'ess 
P2s-sinperi*cosascnd-cosperi'*sinasEnd' cos 
P4=cosperi*sinascnd+tsinperi*cosascnd*costilt; 
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Sin Der sinascrna+cosperi*cosasend*costilr; 
El Si ee ean Samt il t ; 
Er =ceosper:xsjınLilt; 
PELLI EE EE EE OR RK UA ARA AAA AAA xx 
> A AE des 
IA AREA U I I U FR TI RK U TE TI TI TU I TU RT RK KK KKK K 
oneJan92=98.920250931*pi/180; 
TUperday=107.0879334097483; 
ao pPerTU=806.81359; 
ME... npels (3); 
time=inpels (4); 
Miao —Tix(time*le-4) ; 
Man=-fix(rem(time, le4) *le-2); 
sec=rem (time, le2); 
nourer-2ihrs*3536004m2zn^*60-sec)/secperTU; 
Bonvch-rem(earchrot*Tlperday* (Tfilter/TUperday+t... 
m0 cneJjan97,7*pa2y/earthrot; 
NH d-T-ilter; 


Tkeep (1)=JD+Tfilter*secperTU/86400; % Save oriciral time 


p 


EDI RE ONU ON OK OR RIO ON OE KAR AR EE U U N U FU U I I TI IT IR TU TR TI TU N 


Ti tee Ze N ike OF FIRST OBSERVATION 1 


DEDI I UU XN OK OK EDE EER ERA or Ok OSEE ET EER REE EE TER R REK 


time=observ (1,2); 
Ius-rix(time*le-4); 
min=fix(rem(time, le4) *le-2); 
sec=rem(time, le2); 
pueobs-(hrs*3600-min*o0-sec)/secperTU; 
Sk X OX CK OO OK OK OO ck Ok ok OO OX OCC OK IA ARA AAA AAA RAMA A AA UX 
n oe N EN OE ENIGE BANE COORDINATS OF SATELLITE POSITION * 
RE KRAKE OR KA RE EK RR ASE RR ERROR REK RE REK RE ER EER EER KEER EER KERE RR K 
theta-twopi(earthrot*Tgnwch-lonx); 
sintheta-sin(theta); 
costheta=cos (theta); 
meence=.00314+(-sinfinc*costheta’*Pl... 
=sinfinc*sintheta*P4... 
NE E COS TADA É Caleulate current magnitude 
ae meine gescherar’r2...2 ef out-of-plane coordinate. 
EMM ie sintheta PS 
Doce WC*FPB)*r*e6inTA: 


SE Ck ck Ck Ok OK Ck Ce Ck Ck Ck Ok Ok Uk C ke e ck ke e Ok OK oko Ok Ok OK OK OK Ok c A AA A AA A AAA A AA AA A A A 5 


EE Oo: GAT ORBITAL ELEMENTS TO NEXT FENCE-PLANE INTERSECTION * 


So RA AAA AAA AAA A A ACA A A A AAA A AAA AAA A 


EZ—1.082645e-3; 
pl, 
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-0; 
for K=1:K 

I KOoop-pr/erzmeown ME. 

GK KKK RK KK KKK KKK KKK KR KKK KR KEK RR KK KK KK RRR KKK KKK RK KK X 
= INCREMENT MEAN ANOMALY & APPLY SECULAR N ie) x 

Gow KKK IIA A A A A A A RR ER RE ER ER ER RR RR RR RR X 
mnmotion xon I E 

xnow (6) =twopi (xnow (6) +mnmotion-Tleor 

EA (k+1l) =enewton (zZnor (6) ,xnon (6) ‚snov 0 


cosTA=(cos (EA (k+t1) ) -xnow(2))/ (1-znow(2) "cos (F2Ar ua ee 
sinTA-sart (1-xnow (2) ”2)*sin (EA(kt1)) / (1-xnow (2) *cos (EA (KEER 
rexnow(1)*(1-xnow(2)^2)7 (l*xmow 2G cM E 
A(k*1)-twopi(asin(sinTA)); 
if TA(k+1)<pi 
ih oO SrA 
TAK LIED RES % Perform dguacrans es 
end E for true anomels 
else 
if cosTA<O 
da (EYES TE E 
ena 
enc 
Br KKK KKK KK RRR KR KK KR RRR KK RK KKK KKK KR KKK KKK KKK RK WK KKK KK KK KK KKK 
2 APPLY SECULAR VARIATIONS TO omeca Sc OMEGA ^ 


GARA AAA AAA AAA RIA TI CK IK TI TI TI TU HK WK KR U U I U KH RE 
delTAsec-turori (TA (ki TAC): 
semilat- “nov (il) snow AI 

Old sem laco. 

SOMEGA- So oe 

Somega-(-01)*(2-572*Ss10N 99 300 200 E 

A SOMEGA*delTAsec; 

delsomega=Somega* E 

xnow(4)-xnow(4)-*delsOMEG 

«now (5) =xnow (5) +delsomega; 


GARA CA A AAA A A A A KK KK KK RK KKK KK KKK KKK KK KKK 


* CALCULATE ATMOSPHERICUIDEUS MD 
GRKK KKK RK KK KKK RK KKK KKK KR KKK RR KKK KK KKK KK KK KR KK KK KK KKK RK KK 
k2=3e-6; 
RS cos: 
atmdaens=k2*exp((l-r)/k3); 


GARRA AAA ARDE A A % 


x APPLY DRAG PERTURBATIONS TO a & e =* 
GK KK KIA A AA A A A AA A A A A A A 
k1=-4*6378.135; 
deldrag=[k1*xnow (1) *2* (1+xnow (2) ) “atmaens ae rise. 
k1*xnow(1)*(1-xnow(2)^2)*atmdens*delTAsec; 
oor ODE 


snow=:now+delarag; 


REAR GEE ULL UST EED EE ER EE RE EER EER EE ER 2.2 22 2 2 2 2.5. 
E CI SSESBBESTRANSPFORPDATTITONOMATRYX (ORBITAL TO INERTIAL) 5 

TR N UILLAM eee EERE Ee RK eK KK 

cosperi=cos (xnow(5) ) 

sinperi=sin(xnow(5)) 

cosasena-cos (Xnomw (4) 

sınasend=-sın(znow(3) 

) 

) 


r 


% Calculate sines & cosines of 
$ angular orbital elements. 


o Ve 


7 
7 
) 
) 
sos - cos(zmew(3)); 
sintiltssin(xnow(3)); 
Pl-cosperi*cosascnd-sinperi*sinascnd*costilt; 
P2=-sinperi*cosascnd-cosperi*sinascnd*costilt,; 
P4=cosperi*sinascnd+sinperi*cosascnd*costilt; 
po--osyoerasPmudseausssspendseosascnd*costilt; 
P7=sinperi*sintilt; 
P8=cosperi*sintilt; 


MT een nur RE EE PORK RR RR DRR 


Š UPDATE TIME = 
EE EE ke Kk ek ek eK KKK RRR KR KR RARER RE KKK KKK KKK KK 
Dann eceh-Ignwen.Tleop; 
Rei liver-If)]lter.?7100p; 
RR soa RI KR ee eK a 
ms ET ME GO TOE EENGE -EIANE COORDINATE OF SATELLITE POSITION * 
ZAR IIA IIA AA AAA AAA AAA AAA AAA AAA AAA 
RR Ma IÉpwopa(earthrot*Tgnwch-*lonx); 
Sintheta=sin (theta (k+1)); 
costheta-cos (theta (k+1)); 
zfold=zfence; 
arence=.D0s1+(=sintine*costheta*Pl... 
SS Mie” Sinthete* Pdo. 
ASS SEE COSTA oi $ Calculate current magnitude 
NS. es net az oz out-of-plane coordinate. 
zsınsıme-sinthera’Pp5... 
FEOS E ES) r SinTAh: 


E IR TE HT TH Sk slc kc eo ok ok oc oc cc ck X Sk ok ok oe coc ok occ ok ok Ok OX OX SK Ox CK OK OO Ok OO ok ke ke e 


NE c D FOL EENEE PLANE INTERSECTION -* 


SAHARA AAA A AA AA A A A AA A A A A A A A 


ME NEGATIVE TO POSITIVE ?? S 
SS ck AA A AA AA e Oe oe oc ok A A A A kc ok RR RR RR RARR RR RR RR RR RR RR RR OE RR RR RR RR ok 
if zfence>0 
uA cs 
n=n+1; 
negtopos=n; 
xnow-xnow!; 
mie med Dier Tonweh, Sao ow, 2, COSTA, ..% 
SinTA, TA(k+1),theta(k+1),EA(k+1),zfenceJ=... 


83 


inttemp(xmo OSTA er 
SinTA,TA(k*l),theta(k4*1),ER (XO) MER N 
P5,P7,P8ê,Tanwen, Tiilter, TIK 
KAOWSKAOW E 
end 
end 


KKK KK KR KK KKK RK KK KK OK OR KK KK RK KK KK OK KR AA AAA A A A 


* POSITIVE TO NEE MITE HAN: 

Ss Ck ck Ck Ck ck OK Ok OO Ok AR RR KOR REK KERR RAKA REKE EK RR ERROR RR RR RR RR RR RR ER RR RE ER AR RR AR 
if zfence<0 
E o ad 

n=n+1; 

Postoreoa 1. 

xnow=xnow'; 

[Tloop, Tfilter, Tanwch, Flag, snow, r eo eee 
sinTA, TA(k+1) theta(k+1) EA(k+1) zien] MA 
Iinktene(znow,r, cover 
sinTA,TA(k+l) ,‚theta (k+t1) ,EBA(k1) BA ES TAN 
FP5,PI,PB5,Tonwch,Tfilter qi oss 

xnow-xnow'; 

end 
enda 


TAOLG=TA(k+1)2 
BRAK KKK KKK KR RRR RR KKK KKK KR KKK KK KKK KKK KKK KKK KKK RR KKK KK RK KR KK KK KKK 
A SATELLIITE HY NAVSELSWE RE Ns = 


GRR KK KEK RK KKK RRR KKK KKK RK KKK RK KKK KKK KK KK KKK CF KKK EK KR KKK KK KKK KK KKK 


if flag== 


SRK KK KKK KKK KK KR KKK KKK KR KKK KKK KKK RK KR KKK KKK KKK KKK KKK KKK KKK KK KKK KK 


E CALCULATE. TIME FROM LAST FILTER UPDATE (ASM GE) = 


BRK KER KK RR KERR RK KKK KKK KKK KKK KAR TI TH KK TH KK I TI I TH KK I TI I TU TI I KK OO OO OO 


delobs-timeobs-Tfilter; 
while abs(delobs)«1 
ifctlength (obser Cr ULT 
delobs-timeobs-Tfilter; 


delT=delobs; 
delTint (p)sdelT*secperTU; 
Tint (p)=(timeobs*secperTU) /3600; 


Tupdate=timeobs-Told; 
if Tupdate<0 
Tupdate=Tupdate+TUperday; 
ena 
Tola= times: 
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Ta sT htertdelT: 
Tgnwch=Tgnwch+tdelT; 


EN U RI I WU IR EI U FI KK U U U TU I I U U I U N N A U FU k 


BEI NCEOF ESF ETF ZEERIODIE VARIATIONS E 


S kk Ok Ox Ok Ok Ok Ok Ok Ok OO Oe kk kc kc Ok RRA AAA AAA AAA AAA AAA RR A A KA 


xmedgr-2mow,. 


daper MA “Er (Ii) /-(2~sema Nae Ze (t=xnow (2) °2))*. 
So Sees cos A mow to) +TA(kK+1))) > 

deper co AE EEN 2 Ms cemoxnow(3)^999weosTA-.. 
M im ue Ses" mow UN TACKCHI)) t. 
As Img m 0) 2*c0s0 7 snow (5)4+5" TE (k+1)) ) > 


daper= 32, 2/12: semrlaerr2)*sin(Z’rnow(3))in.. 
es ao TAKEL) 
GeOMser=37 07) (Aone ae Zi cos(xnow(3))*... 
sn (ae TA ELI) 
en ea WA eme 2) IE ET SIN (HAD (ADILY se 
Sn Ma o a TA (LIO 
ee me- AE (2 sem flar 2) = (1/xnow(2)*. 
aux 2 ears (AKSIE TAS a. 
is Ed EG sin (2 *ROWIES TART ee 
oe Gr ae (N 22 snn(2°2mow (5)=+3*°TAlkrl)))r... 
mr canon) 2)*sàân(2*TA(k+1))- 
(1-5/2+sin (xnow (3) ) 12) * sin (2º (snow (5) VTA 11) ) 
ES Rs mo nor 2*e.srn(2*xnmnow(5)44*TA(k*1)))); 
comper=dtemp; 
dMAper=dnewper-domper; 
varper-[daper;deper;diper;dOMper;domper;dMAper]; 
xnow-xnow-varper; xnow(6)-twopi(xnow(60)); 
EA(k-1)-2enewton(xnow(6),xnow(6),xnow(2)); 
ai (ees (EA (EI) Snow (2))/(L1-=:now(2)*Cc0s (E 
EEN set (1 mor (22) *sIDA(EA (kT1) ) / (1-2now 
Een dir (AE / (1+xnow (2) *cosTA); 


mu 2 
EI ee GE LS 


EK I)etwopitassm(sznTA)); 
AA pL 


iN EST EG 
PL EA OE: $ Perform quadrant check 
end % for true anomaly. 
else 


if cosTA<0 
TAGS N pPI TAk]; 
end 
end 


Bo KR KK RI Ik A A AA AA AAA A AA A A AA A A A 


one o> tee ANGULAR DISPLACEMENT OF RECEIVER SITE FROM I-AXIS * 


Bok ARA AA 
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phi=twopi (earthrot*Tonwch+1lon (observ( pa eae 
sıinphl-siınlenss 


sinlat=sin (lat (observ( 
coslat=cos (lat (observ ( 
sinlonssinileniess 
coslen-cos (long er ses, 
cosperi-cos (now 057757 
sinperi-sin(xnow(5)); 
cosascnd=cos (xnow (4) 
Sinascnd=sin (xnow (4) 
Fr 
Er 


U'U'U'O 


Costrlt- Cosimo ME) 
Sint i lt—=sint RO IS) 


7 


; 
Pi=cosperi*cosascnd-sinperi“sinasend sos ERT 

P2=-sinperi*cosascnd-cosperi"sainasence? cost 

Pá=cosperi*sinascnd+sinperi”cosascne cost 

P5=-sinperi*sinascnd+cosperi*cosasena- cece mar 
Fiszsinperi' smelt; 

P5s*ecospers SINNE TYE: 


sinaz=sin(az (observ(p, 3))); 
cosaz-cos (az (observ (>, 3) )- 


Hiscosias cos: 

H22cosTtat"' sym 

H3=sinlat; 
H4=-sinphi*cosaz-sınlat *cosph72s 7 u 
Ho =Ccosaz*cosphi=sinaz*sintat sin pe 
HE-simas cos llar 
A7=sinaz*sinphi-=cosaz*sinlat*cospuis 
N8=-sinaz*cosphi -cosaz*sinlat?sıy me 
H9=-cosaz2eosstat: 


R1sHa Phir Pa He BY; 
h2-H4^Pp2THO5 EE LUG DE 
RAZHJ FREI TSHYZPAL eas 
RS-HI/SEZENSB FSLN, 
RI HIEI MZ Dn 
RöS-HIZPZEHN 2 PS Ea 


LARA AAA AA A AA AA AAA AAA AAA A AAA AAA AAA 
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* ESTIMATE MEASUREMENT 2 
A A TE KK IK TU TR TI TI KT U U TU I KK TU I IF IH KH N U IT TU KK U N X 
Bee m29057200729607720*si1inlat^2)-*alt(observ(p,3)); 
röffset=.003]; 
rhoH=r* (R7*cosTA+R8*sinTA) -Rsite; 


THOB- FA (Riese oes TAYR2*siınT2); 

rhoN=r* (R4*cosTA+R5*sinTA) +roffset; 

so Saa non 22: DoE arnon 2); 
Reco (Dj (rhokß/rno,chel,rno]; 


size=xnow (1); shape=xnow(2);tilt=xnow(3); 
ascnd=xnow (4) ;peri=xnow (5) ;MA=xnow (6); 


PEEL | nn AR EL ek nes hn Ke KKK KR RRR RK RK KR kK k 


= CALCULATE COVARPANCENMATRIM P x 
EE to E Ii OE AE OE ARE AE AE RE KERR AE KERK RR RR RR RR RR NN 
J2=1.082645e-3; 
Sins A=sintEA (kt la 
SOseh—-Cos (2A (k+)) 7; 
sindecl=sintilt*sin(TA (x+1)+per1); 
semilatzsize*(1-shape^2); 


LARA AAA AA AAA AAA A A A A OR RR CO CAO CAO A OK OK OK OX Ck Co 


* CALCULATE APPROXIMATE delta(TRUE ANOMALY) É 
ERRAR AAA RA AAA DAR RK RK KR Rk kk ao kk Rk RK 
TOE OO now (])” (AA); 
delMA-mnmotion*Tupdate; 
delTA=delMA; 
ls eae 2 SiNEAT LT shape SINBATCOSBA; 
A2phi-1-shape*cosEA; 
Dx as 2-S AA sint11E RI: 
Zaehnı-1-3*sindecl”2; 
Loonsz-ygj-shape^2; 
A6phy-sdgrt(l/size*J2/r^3*A4phi); 
A7phi=1+shape*cosTA; 


dOMEGAda-3*J2/(size*semilat^2)*costilt*delTA; 
dOMEGAde--6*J2*size*shape/semilat^3*costilt*delTA; 
dOMEGAdi-3*J2/(2*semilat^2)*sintilt*delTA; 

EOoEEengM-coI» M 2*cemilat^2)*costilt*AlIbohi/ (SinTA*A2pni^3); 


domcgada PS ZE “semi lat 2) *A3pni*delTA; 
domegade-6*J2*size*shape/semilat^3*A3phi*delTA; 
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domegadi--15*J2/ (2*semilat”2) *sinti Me Eosti MM 
domegadMs-3*J2/ (2*semilat”2) *A#phi*Alphi/ (sins EIN 


drkdak=A5phi/A7phi; 
dMAda=3/2*Tupdate*A6Gphi* (-1/size”2+J02*A4pon1* (-37 2 eee EN 
drkdeks (-2*size*shape?* (1-shape*cosTA) -semilat*cosTA) /AIpm 
dMAde=3/2*Tupdate*A6phnhi*J2*A4phi* (= aards 
dsindkdi=sin(TA(k+1)+peri)*costilt; 
dMAdis3*Tupdate*Afphi* (-3*J2/r”3) *sindesl as MEE MEE 
dsindkdo=sintilt*cos (TA (k+)) per Af 
dMAdomegas3*Tupdate*Afphi* (-3*J2/r” 3) *s1indecl derd cE 
drdTA-semilat*shape*sinTA/A7phi^2; 
dsinddTA-dsindkdo; 
dmessdTA=-3*J2/r*2*drcTA- (-9*J2/r 4" sindec Zac se. 
6*+J2/rºs*singeclasincan 
GTAdMA--Alphi/(sinTA*A2phi^3); 
AMAdM-143/2*Tupdate*Afphi*dmessdTA*ATAAMA; 


dada=1+2*ki*size™ (l+shape) “aumdens mar 
dade=kl*size”2*atmdens*delTA; 
dadM-z-kl*size^2*(l1*shape)*atmdens*dTAGOMAS 


deda=ki*A5phi*atmdens*delTA; 
dede-1-2*k1*si2e*shaperacmdensrgert 
dedM--kl1*size*A5phi*atmdens*GTAdMA; 


phimat-[dada dade 0 0 0 dacM; 
deda dede 0 0 O dedM; 
OO 21070 020% 
dOMEGAdGa dOMEGAde dOMEGAdi 1 0 dOMEGAdN; 
domegada domegade domegadi O 1 domegadM; 
dMAda dMAde dMAdi O dMAdomega GMAGdM]; 


KA RARR RARR RR RR RR RRR RR RR KRRR RR RR RR RR RR RR ARE RR RR RR RR EER RR K RR EE EEN ON 


% 

x CALCULATE PLANTINOISE x 

SEC Ck OK Ok C CO OO Ok Ok A AA AA AAA A A A RA AAA AAA A A A AA OR KERE EER RE RR RR RA A AA 
J2=1.0826e-3; 
semilat=xnow(1)* (1-xnow(2) 72); 


siga2-(3*J2*xnow(1)/(2*semilat 2*0] ome a ee 
sin(xnow(3))^2*co058(2*(xnow(5)s Th RE 
sige2=(3*]2/(2*semilat“2)*((1-3/2*sin(xnow(5))' 2) aa AA 


s1gi2s (3*J2/ (8*semilat2)*sin (2 ow 
cos (2* (xnow (5) TAK I 
mnmotion=xnow (1) * (-3/2); 
MApert=mnmotion*3*)2/(2*semilat”2)*(1-3/2*sin(Xxnow(5)) DE 
sart (1e 2 
SigOMEG2- (3*J2/(2*semilat^2)*mnmotion*cos(xnow(3))*Tupdate 7 
sigomeg2= (3*J2/ (2*semi lat 2) mamona na 
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MISS MAS TAO LS) 
S:gMA2-2 (MApert*Tupdate) 
e- [eS 2esigaz 0-0 0 
De sides 00 0 
Sgt Ze Cn 
0 sigOMEG2 0 O 
0 0 sigomeg2 0 
OTORO SIOMAI Te >> 


op EE date) 2; 


/ 


2 
0 0 
0 


CO OC Ca GE 
OO EE 


Ms OO OE OX DE RE KAR KEE ERR KKK KKK KK KK 


= PROPAGATE ERROR COVARIANCE  * 
MEE a XA ELE KDE X A A ee EO EE KERKE EA RR RR KI KYK RR IR RR RR A 
Pkk=Pkkml; 
Pkkmi-phimat Pkk phimat "tO; 


ER ER I RU I U I N I I I KH U IE I KU I I I RA E RR A RÃ k 


n CA Nn FOBSERVATION DATA = 
EE EER IH I IE RT TR A IF e A or X KK o o N 
cosalpha=sin (observ (p, 4) *pi/180); 
cosbera=sin (observip, >) *pb1/180) 
zZ(:,P)=(cosalpha; cosbeta] ; 
per NR ee ONDE I EL I NUN I N U N TI KL TH AE KEER OE EO KA U I TE TI IK U RR I I U I U I E KOR ROK ak 
E CALCULATE LINEARIZED MEASUREMENT MATRIX H E 
EA KRAAK EE DD RADAR RA O NT TU NK TE PT RN CK N TH IT HT RO RKO RK KK RK AK A 
al=]=e2 non (22772: 
az = 1 xmow (2) CoS TA: 
a3=2*xnow (2) +CoSTA” (1+xnow(2)%2);> 
at=2 “know (2) SINEA*~coshA; 
dU 2 no ZS 
a6=-a5*sinEAta4; 
a7=1-xnow (2) *coSEA; 


erda-al/a2: 
dere Il “aap: 


dcosTAdT--sinTA: 

dTAGEA- (-a5*sinEA-*a4)/(-sinTA*a7^2); 
dEAGdGMA-1/a7; 
dcosTAdM-dcosTAdT*dTAdEA*dEAdMA; 


dsinTAdT-COSTA: 
dsinTAdM-dsinTAdT*dTAdEA*dEAGMA; 


educ mou al*SintA/ac 2: 
drdMA-drdTA*dGTAGEA*dEAGMA; 


elle sinperi'sinasendisinti ERHAS 


Sinperi cosa Sond i Sintilt*H5+.:. 
sinperi*costilt*H6; 
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" 


aRZ2di=cosperi*sinascnd=sintil 2 ma 
cosperi*cosascnd^sSint2 ND OM 
Cosper m cosie HE 
dR4di=sinperi*sinascnd*sintilc mie 
sinperi*cosascnd”sinerte er NM 
sinperi costi li N 
ARSdiscosperi*sinascnd'sinti le MD 
cCosperi*cosascnd' sinti ma 
cOsperiicost VET HE 
dR7dis=ssinperi*sinascnd“sint ms 
sinperi*cosasend’sinerm.e 
Sinperi costs HI 
ARBAdi-cosperi"*sinasend” sine eS M 
cosperi*cosasend”sınc ı Le 22 a. 
Cösperi’costj ein 


GR IAOM=--PA= A PIE 
ORZdOMcSePDESHMA LE MS 
ARAGOM--EBA HI ED SEE 
GRSAOMS-B EP RAY EP 
GRIAdDOME-PA acl: 
ARBAOMSSESHI DE EE 


dRldomzP2*H44P5*H54PB8*H6; 
dR2dom--Pl1*^*H4-PA4*H5-P7*H6; 
GR4dom=P2*H7+P5*HB+P8*HSG; 
dR5dom--PI*H7-PA*HB-P7*HO9- 
arR7dom=P2*H1+P5*H2+P8*H3; 
aekbdom- - PI=HI=-PAFTIZ- E72 1: 


doHda-(Hi*cosTASRBSO^S3DnDA) vods. 
ODEEas (RI '*EOSTA Rosina ce) 
apNda= (R4*COSTAFTRI*SINTA da: 
dpda= (rhoH*dpHda+rhoB*dpEda+rhon “Genoa wae. 


dbpHde-(R/*cOSTA-RNS*SqnrA)jemde 
APEGes (Kl *COSTAR2 SIR sies. 
dpNde- (RÁ*COSTA-R5*sinTA)*drde; 
dpde= (rhoH*dpHde+rhoE*dpEde+rhoN*dpNde) /rho; 


dpHdi=(dR/di*cosTATGR@G@i~< iia 
GpEGi=(dRidi*cosTA+eRZGqi sip 
ApNdi- (dR4Adi*cosTAARSA1 *s ir Ba) Er 
dpdi- (rhoH*dpHdi4rhoE *dpEdirhoN asma mae, 


APHAOM- (dR7dOM*cosTAdRBAOM*sinTA) *r; 
APEdOM- (AdR1dOM*cosTAAR2:@M SIn 


90 


dpNdOM- (dRAdOM*cosTA*dR5dOM*sinTA)*r; 
el ele (EER deHAOMrAsE SdEEAOMRON*dPNdOM) /?Zho; 


derem (AE ASM SOSMA ARSAOM SINTA) *E; 
ape dom (dridomícosTAraRZaomsiniTA)*r; 
dpNdom= (dR4dom*cosTA+dR5dom*sinTA) *r; 
gdecen— (sok doHASmtEROE dEEASMI HON dPpNdom) /rho; 


eke dcos BE AsindAaM) *r+(R/*cosTA+R6* 
SINTA) draMA; 
Genova RA dos HE Ba GE nTAdM)] r+ (R1*cosTA+R2* 
sim UA adn amas 
dpNdMA= (R4*dcosTAdM+R5*dsinTAdM) *r+ (R4*cosTA+R5S* 
sinTA) *drdMA; 
dpdMA= (rhoH*dpHdMA+rhoE*dpEGMA+rnhoN*dpNdMA) /rho; 


drneokeao  rEde/rhe rhob*eapda/rno72: 
cueldasdesdossuserhoNsdpda/rshHo^2; 


di gomde-desbdeiuscnhoEssdpde/rhe 2; 
drhoNde-dpNGe/rho-rhoN*dpde/rho^2; 


eiersel se io HOE dpdi/hof2: 
Sunclur-deNdi)rhocorüoN5dbdi/rho^22; 


duuns5sdgdoOMzdbpsEdoMnhoerhob*dbdOM/rho 2; 
dari NdoOzadpNuoeMAmbo-rnoN*dpdOM/rno^2; 


dnhosdomsebEdom/rho-rhoE*dpdom/rho^2; 
drhoNdom=dpNdom/rho-rhoN*dpdom/rho%2; 


OpionaMA-cpnoMA/  rho-rhok ~cdpaMay rho” 2; 
drhoNdMA=dpNdMA/rho-rhoN*dpdMA/rho%2; 


D suesd-EcdrHobEde drhoEdri'drhoEdOM drhoEdom drhoEdMA; 
drhoNda drhoNde drhoNdi drhoNdOM drhoNdom drhoNdMA]; 


S5 Ck CO Ok C Ok C OC OO Oc Ok CO CC COO OO CO OO OK RÃ AAA RA A 


* CALCULATE KALMAN GAIN * 


Sc ck Ck Ok CO OCC Ok OO OO Oe OO Oc kc Ok Ok Ok Oo Ok Ok c Ok Ok Ok Ok Ok Ok III AR ER RE RR RR EER RE RR EE RARR RR ok 


GCSPKKMI H *inv(R*Pkkmi*H'+ŤR); 
Restelo) Zi pl-zest(t,P); 


S. Ck ck ck Ck C AAA AAA A K 


“N UCRTE IMPROVED COVARIANCE * 


Sc ck Ck ck Ck Ok Ok Ok Ok OO C OK COO Ok Ok Ok 0k ok Ok ck Ok Ok 0 ok Ok Ok Ok Ok Ok Ok Ok Ok Ok ck ck ck Ok Ok Ok C OO Ck Ok OX Ok Ok Ck Ok Ok kk ko ox 


Pkk= (eye (6) -G*H) *Pkkml; 
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E * 


€ 


a 


PkeepD (:,D” 65 le plc 


* RR RR ERK RR AR RR RR AAK RR RE R RE RR RR REK RR ER ER RE RR RR RE EE ER EE RR ER EE RE RR E X X 


CALCULATE IMPROVED MEAN ORBITA TECE T = 

X* OK ok OK CX OK Ok OK CK CK CX OX Ok OK OK OX CK OX CK Ok CK CK CK OX CX OK CX CK Ok Ok OK CE C CX 0 OO CK Oc CK OX CK CK OK OK OE OX CK CK CK CK CK CE CX RR ERK RE OX CX OX OX 
xnow=xmean; 
Gresijei zer 
xnow=xnow+G*resid(:,P) 
xnoOw (6) =Ewore (ae. (ee 


we "oe 


xX3int llo 2 Ao DO 
KXint (3367p) =xnoOw (so) 3 6 a, 


C. * x X* x ook ok x c ok ck OX OK OK OX OK OX OK CK ok Ok OX Ok C CX 00k Ok OE TI KU TI TO I TI TI UOTE HOT CK OE CX CX CX OX 0X 


CALCULATE IMPROVED ECCENTRIC ANOMALY * 


RR RR RR RAAK RE RR RARR RR RR RR RR RR RR RR RR RR RE RR AR RR RR RE RR RE RR RE EER ER EER ERK RR 


A(k+1)=enewton (xnow(6), XnNOW(6) , xnOw( 2) 


* % RR AR AR RR RR AR RR R RARR RARR RARR RR RR RE RR TERE RR RR RR RR RR RR ER A E E 


= SOLVE FOR IMPROVED cos(TA), Sin(TA), z, Si idea at 
RR RE AAA AAA HATE xx x xx 
Son A(k+1)) -xnow(2))/ (1-xnow (2) "eos (Ei (ka 
sinTA=sart (1-xnow (2) °2) "sin (EA (k+ 1 Al A 
cose ER G1); 
PS Noli Nel 


Jr coc 


AXk-l) oi.r2 2 + Perform quadrant check 
and * for true anomaly. 
else 


if cosTA<0 
Taki) = +p TA 
end 
end 


TAOIG=TAIk I 


SE Ck OK OK Oc kk ke eO RR RR RR RR RR RR RE RR AR RR RR RR RR KERR RR ER EER RE ARE EE A 


O osea Mm 
p-p*l; 
time=observ (p, 2); 
hrs-fix(time*1le-4); 
min-fix(rem(time,le4)*1e-2); 
sec=rem(time,le2); 
timeobs=(hrs*3600+min*60+sec) /SecperTU; 
delobs-timeobs-Tfilter; 
else 
delobs=2; 
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end 
end 
end 
end 
flag=0; 


S.A ER EO EED AE EI HK ER KAAK RR EE REK KYK RR N x 


fees on! TIMES POR NEW DAY > 
Ex Xxx A eee KARR KK AEA ERA KKK KE KR 
Ar “arwER 10 1E EI ES 
Reanwch- lenwere lo SS Sos, 
end 
Er NON OS cee 
sta ie ie OE el iese OE Oe IE y 
JD=JD+1; 
end 
xkeer (:,k+1)=xno0w; 
Tkeep (k+1)=JD+Tfilter*secperTU/86400; 
end 


In-plane Iteration 


So ere olie, TanwEeRk, Lag, snow, r, COSTA, ... 
Er LI” cnera, Bnhow,zftence|]= as 
Tor ven ne cosTA, sınTa, TA, theta,EnsW, :.. 
EIE EE EA Ee TEnwEh, Tillter; TT): 
DEOW-ZDOWw':; 
J2=1.082645e-3; 
LEorcthrort-.05883378]771654;*tDefine earth rotation rate(ragd/TU). 
Mon :=-"101.31348*p1/180; 
MA e-— 5 .58310*pi/180: 
ERR t inc-sain(finc); 
Basftine=-coös(finc); 
delT=0; 
zfence-1; t$ INITIALIZE OUT-OF-FENCE-PLANE VALUE OF POSITION 


CARALHO A A 


% iS EIS TE TOU BLANE CONDITION "= 


GRAAL A 


while abs(zfence)»-1e-8 
SE Ck C A E I KK TH A A A A A A AA A A A KK KK KK KKK KKK 


% = ADD IN PERIODIC VARIATIONS n 
Swe wk RK KR RARA AAA A 
semilatexnow(1)*(1-xnow(2)^2); 
dor. ME ie ON (Oo sem aE IE (IE KAOU II AI IE 
Sino ers Cos (Z2* (xnow (5) +TA)) > 
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deper=3*J2/(2*semilat%*2) * ((1-3/2*Ssin (know (3))) - 2) M 
1/4*sin(xnow(3)) *2*cos (2Z2xnow Cn 
T/1l2*sin (znow (3) ) *2*cos (Zuger eee 
diper=-3*J2/ (8*semilat”2Z) *sin (22 nem ss eee 
"cos iz zuew(S) + TADE 
dOMpers3*J2/ (4*semilat”?2) *cos (snow (3)) sin (2% (ars RE 
dnewper-s-3*J2/ (4*semilat?2)* (1-5/2*sin (nor N 
* Sim (ZA (now o M E 
dtempsz3*J2/(2*semxlat^ 2) (1 /xmnow0 DM 
((1=-3/2*s1in (snow (3) ) 2) “cando 
1/4*sin (xnow (3) )°2*sin(2*xXnew (Ss MD 
7/12*sin (xnow (3) ) °2*sin(2* know (5S) N 
1/2* ((1-3/2*sin (xnow (3) ) 7 23 AS (Ze M 
(1-5/2*sin(znow(3))?2)*sin (2 (nor (ar EN 
3/8*sin(xnow(3))^2*sain(2*xnow(5) SINO 
domper=dtemp; 
dMAper-dnewper-domper; 
varper-[daper;deper;diper;dOMper;domper;dMAper]; 
xper-xnow-*varper; 


EAper=enewton(xper (6) ,xper (6) ,xper (2)); 

cCosTApers(cos (EAper)-zper(2Z))/ (1-sper (2) *sosia E 

sinTAperssart(l-xper(2)”?2) *sin (EAper) / (i-xper (aa 
“EOS (Ape 

rper=xper (1) * (1=-=per(2)72)/(1+=per (2 costa ae 


TAper=twopi (asin(sinTAper)); 
1 TJAper<pi 
if cosTAper<0 
TAper=p1-TAper; % Perform quaGvane pe vee. 
end % for true anomaly. 
else 
if cosTAper<0 
Täper=-3’np2 Tier 
end 
end 


E CACA CA A CCA A AA A A TI TI TI TI TU A TI TI TI I TI TI TI TI TI TI TI A A AAA A AAA CA 


% 
% x RECALCULATE PERIODIC TRANSFORMATION MATRIX 2 
ZARIAAKIAA KAKI KI KK KAKI AAA AA AAA A AAA 
Piper=c0os (xper (3))*cos1xper (AMAN 
sin(xper(5))*sin(xper(4))*cos pem ll 
P2per=-sin(xper (5)) 26088 (spez 2 
cos (xper (5) ) *sin(xper (A)) "eos Sperre ide 
P4per=c0os (xper(5)) *sin (zper nn c 
sin(xper (5)) *cos (xper (4))*co0s (perrea 
P5per=~sin(xper (S) ) *sinixpen Tm 
+cos (xper (5) ) *cos {xper (A) cos (Pedi E 
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Biper em peeis))*sin(xper(3)): 
Föper=-eos(xper(5))*sin(xper (3)); 
Sure em mereartchrot*sper (|) l.5*rper*... 
cosTAper* (Plper*sin(theta). 
mo pe scene jcheta)jscinfinc* eartnrot*sper (1)°1. oS 
mee sin EST (PPPETTS N (Eheta) -P5per”*. 
cos (theta) )*xper(1)/(1-xper(2)*cos(EAper))* 
(cin Papen ia (ac lM SPUDer*cosítheta)-.,.. 
Sint rne*PAper cim@ (elec aie ecostinc~P?per) + 
sad Era eccle e-sanfine^P2per*.. 
Seoslevera)-S hres: seen cm oEnera)+cosfinc*P8per) ); 


Amer OU N. si EG eos neta) ~“Piper-sinfinc*.. 
cnmdtehetc cL coSPIDO*P7Der)*rper*.. 
ees TEEL (SI ANE CoS theta) *P2Der... 
= romeu checa EDE. EOSLinNC PBpEr)?T 
Poem" sind Ape, 
delMA=-zfence/stuff; 
mumoctconcemowtl)^t-57/52) 
del2T-delMA/mnmotion; 
EX UNO XK A AK KKK AK AA KK KK RR AK RA RAK RAK KKK RK RKRKRKKKE ARERR KKK KER KE KK KX 
Oe Steere E TT AL ELEMENTS TO ITERATIVE VALUE OF TIME * 
SEEK OC KOC Ok COO Oe eO e ke OI Oe e ke ek KH TH KT KK TH KH A KH TU TU TU AA ERK RAR KKK 
xnow (6) =twopi (xnow (6) +delMA) ; 
Enow=enewton (xnow (6), xnow (6), xnow(2)); 
TAold=TA; 
eo (es: (Ere. ow 7L-umow 9) *cos(Enow)); 
EA sr e (1. me (2) 2) sin (EnSW) / (1snOW (2) si. 
SOS (ENOW) >: 
aar NA (1 new AAI now (2) *cosTA); 


TA=twopi (asin(sinTA) ); 


D LODS 
al EoSTASD 
TA=pi-TA; $ Perform quadrant check 
end % es true anomaly. 
else 


di eon 
IA BYE TA: 
end 
end 
RR RR AAR RR RR RR RR RR AR AR RR RAAR RR RARR RARR RAAR RAAR RR RR RR RR RR RR RR RR K 


% 
% = APPLY SECULAR VARIATIONS TO omega & OMEGA * 
SKK KKK KKK KR KK kk Kk AAA kk ko 
delTA=TA-TAold; 
if delTA>pi/2 
delTA=delTA-2*pi; 
end 
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ir delta gen 
delTA=delTA+2*p1; 
ena 
semilat-xnow(1)* (1-xnow(2)^2); 
Ol=-32I2/sene Tq 2, 


SOMEGA=01*cos (xnow (3)) /2; 
Somega--Ql1*(2-5/2*s3m(xmow (30 MP PE 
delsOMEGA=SOMEGA*delTA; 
delsomega=Somega*delTA; 

xnow (4) =xnow (4) +delsOMEGA; 

xnow (5) =xnow (5) +delsomega; 


SEKER ER ER RR RR RR RR RR RR OR RR RR RR RR RR OE REK R RE RE KERR RE RE EER ER EER ER ER EX 


A CALCULATE ATMOSPHERIC DEN Ix E 
GARA AAA A AAA AA AA AAA AAA AA A A ACA AAA AA % 
k2=3e-6; 
jp esq cu Ee 
atmdgenscek2*'exou lp s 
SRK KKK RK RK KR RK RK RK KK KKK KKK RK KKK RR KK RK RK RK KR RK KK 
ie APPLY DRAG PERTURBATIONS TO a € e x 
Sn RRR KKK KKK KR KKK KK RK KR KKK RK KK RK KK KK KR KKK KK KK aa K RE 
ka ES IE MS; 
celarag-[ki*xnow (1) 27 (1 +xnew ee es 
kl^xnow(l)w(le*xnow(tz)9 7 aa 
Om DS 
xnow=xnow+deldrag; 
GRR KKK KKK RK RAK KKK KK KK KK RK KK KK KK KK KR RR RK KK KR KK KR KR KKK 
* RECALCULATE TRANSFORMATION MATRIX (OREITAI TO TE GRA, x 
GRR A AA A AAA TI TU I N TH TH A A AA AA A RX 
Pl=cos (xnow (5) ) *cos (xnow(4))- 
sıini=»now(5) 1 ’sın(znew(A 
P2="sin(xnow (5) )ºcos nov] 
Gos (now (5): 2 san HOW (4 


cos (no NE 
We 
“cos (ee E E 


) Ee 
) 4) 
) n 
) ) ) 
)*cos (xnow (4)5 
) 4) 
) ) ) 
) jo 


P4=Ccos (*now (5) sin (nek MAE 
sin(xnow (5) dos ne E 
ES==Sinixnow(S)) Sin (non E 
cos (xnow lo "COS nod cos nei 
Pj=sin (xnow(5))*sin (nen cin 


Pb=cos (xnow(5))*simn(xnow (Sr 


* RR AR RA RR AA RAAR RR RARR RARR RAAR RR RAAR AR RR RAAR AR RA RR RARR RARR RARR AAR RA 


% 

% * IMPLEMENT TIME GOEIE TOON S 

LAA A AAA A A A AA AAA AAA AAA A 
delT=delT+del2T; 
theta-twopi (earthrot* (Tgnwch+delT) +lonx) ; 


end 
Bo RO ALLA AA AAA ARA ARA AA AH 
E di UPDATE ACTUAL TIME OF INTER F E A 


SE Ck Ck C Ok Oc Ok Ok ok oe ck Occo oe ck ok ook kc ok ok AAA AAA AAA AA AA A 
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EenwchesTonwclsdejT; 
Mir ltcer=Tiiiter+deltT: 
DTrdelTs 

2now=xnon'; 


MIX X wx. Ww» c» KKK KKK KKK KKK RK KKK KKK KK AK KKK KKK KKK KK 


* CALCULATE FENCE-PLANE COORDINATES > 
AC a A a a I I TI TU KK KK KK KKK KKK KEK KKK HK KK KKK 
Veg atmmeteositneta) *Plteostine*sin (theta) *P4+.. 
Si E ECOS TAS Eos tine” cos (theta) *P2+... 
gesainersın(theta)*PSsrsinrine*P8)*r*sınTı; 
Ms encta)*Pli+rcos(thetaj*P4)rrícosTAÃ... 
+(-sin(theta) *P2+cos (theta) *P5) *r*sinTA; 


BKK KKK KKK KKK KKK KKK KK KKK KK KKK KK KK KKK RK KK KK RK AAA AA AA A A AAA 
npe E Ae inm tore * 


REM CNE WRETRER ANY ROVER MAY BE ABLE TO OBSERVE SATELLITE 
Soon E E A U kK Kk kk 
SE N 202, --20s (al 57/180) 
tlagsl ; 
end 
zÍfence=0; 


oo o? oo 


Eccentric Anomaly Iteration 
function [EA] =Enewton (EA,MA, e) 


ge EN 
RAS 
MA=MA+2*pi; 
end 
end 
errE |. 
while abs (errE)>=1le-10 
errE=EA-e*sin(EA) -MA; 
derEA=(=1/1(1-e*cos(EA)))*erreE: 
DAEWOO po (EA+FaelBA): 
MA-twopi (MA); 
end 
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Angle Between 0 & 2x 


% This function will take any angle (in radians mm 
calculate 
& lts equivalent between zero and 2*pi 


function phistwop iS) 


Jf oss 2 a 
while phi>=2*pi 
pPhi=ph3-2>D27, 
end 
else 
while phi<0 
PERI Ed AE, 
end 
end 
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